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ABSTRACT 



A simple, microprocessor oriented algorithm is developed 
to identify, track and suppress bending mode signals from a 
control system’s rate and position feedback signals using 
adaptive digital notch filters. The algorithm can be used 
to suppress bending modes having center frequencies as close 
as one octave above the control system gain cross-over fre- 
quency without introduction of the excessive phase loss 
associated with conventional lowpass filtering techniques. 

A third order model of the trident missile autopilot pitch 
attitude control loop is contaminated with two dynamic, 
destabilizing bending modes and used as a concept demonstra- 
tion model. The algorithm is demonstrated by stabilizing 
the pitch attitude loop in the presence of two bending modes 
with unknoivn gains, damping, center frequencies and rates 
of change of center frequencies. 
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I. INTRODUCTION 



The purpose of this thesis is to develop a simple, micro- 
processor- oriented scheme that eliminates body bending mode 
signals in a missile control system by adaptive digital notch 
filtering of the bending mode resonant frequencies from the 
rate and attitude gyro feedback signals. Although this thesis 
uses the Trident Fleet Ballistic Missile as a conceptual 
demonstration model, the process is applicable to many control 
systems having unknown and/or time varying mechanical reso- 
nances . 

A. THE BODY BENDING PROBLEM 

Body bending is a composite term which encompasses the 
coupling of all structure -borne vibrations from various 
sources into the missile control system through the mechani- 
cal sensors such as rate gyros and accelerometers. The 
effect of the body bending mode signals in the rate and 
attitude gyro feedback signals is destabilizing since if 
left unattenuated, the missile control system tries to com- 
pensate for the sensed body bending rate through the thrust 
vector control (TVC) system thus further exciting the 
bending mode and the inertial or rigid body rates soon be- 
come insignificant and go uncorrected. 

Some of the sources of body bending are: 

1. The missile acting as a cantilevered beam. 
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2. Hydraulic servo resonance. 

3. Structural bending at the attachment points for 
hydraulic actuators, 

4. Beam-like resonance of sub-assemblies such as 
upper stage motors and payload. 

The body bending problem in missiles is traditionally 
approached in one of three ways : 

1. The mechanical sensors are physically located 
at nulls in the lower frequency bending modes and the 
natural bandwidth of the control system and lowpass filters 
are used to attenuate the higher frequency bending modes. 
Reference [1] contains a detailed description of the analy- 
sis and testing required to implement this approach in the 
control system of the lXSM-65 missile. The popular name, if 
any, of the XSM-65 missile is unknown. 

2. The mechanical sensors are physically located 
at a position in the missile where the deflections due to 
the lower frequency body bending modes are in phase with each 
other. Phase shifting networks are used to ensure that the 
bending mode signals are fed back in such a manner that the 
bending modes are attenuated. Lowpass filters are still 
used to attenuate the higher frequency bending modes. 
Reference [2] contains a general description of how this 
method is implemented in the control system for the Saturn 
Missile . 
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3. Lowpass or bandstop filters are used in the feed- 
back control loops to attenuate the body bending mode fre- 
quencies from the rate and attitude gyro feedback signals. 

A general description of how this method is used in the 
Trident Missile Control System is contained in Ref, [3]. 

In the Trident Missile the above three approaches have 
many undesirable features. Location of the sensors at nulls 
or where the body bending modes are in phase requires a very 
precise knowledge of the bending mode characteristics and 
the freedom in the initial design to locate the sensors at 
the anticipated optimum locations. In a solid fuel missile 
such as the Trident, the location of the sensors is pre- 
determined to be located with the other flight control 
electronics with very little latitude to select a null or 
position where the bending modes are in phase. The second 
major limitation is the lack of precise knowledge of the 
bending mode characteristics. The bending mode resonant 
frequencies, gains, damping and phase are a function of many 
factors including time into the flight, attitude, payload 
configuration, skin temperature, and thrust level. For 
example, the Trident Missile is designed to compensate for 
as much as a twenty-five percent variation in thrust level 
due to slight variations in the solid propellant mixture and 
casting while only a three percent thrust variation is 
acceptable on the liquid fueled Saturn Missile [2] . Without 
a precise knowledge of the parameters controlling the bending 
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mode resonant frequencies, gains, damping and phase, the 
lowpass and bandstop filters must be designed on a worst 
case basis and as a result reduce the control system band- 
width and phase margins much more than if a narrow notch 
filter had been used at each specific bending mode frequency 
to remove the bending mode component from the rate and 
position gyro output signals. 

As the Trident Missile evolves into the next generation, 
the anticipated bending mode frequencies will become much 
lower and approach the natural frequencies of the flight 
control system. The use of conventional lowpass filters 
under these conditions would produce so much phase loss that 
the missile would become unstable. 

B. THE PROPOSED SOLUTION 

The proposed solution to the bending mode problem is an 
extension of preliminary work done by Lockheed Missiles and 
Space Company [3] . 

If the bending mode resonant frequencies can be deter- 
mined in near real time during the flight and tracked as 
their controlling parameters change, then adaptive notch 
digital filters in the rate and position gyro signal feed- 
back paths can be implemented to attenuate a selected number 
of bending mode signals that are contaminating the desired 
rigid body rate and attitude signals. 
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The difficulty in the implementation of such a solution 
lies in the in-flight determination of the bending mode 
frequencies rapidly and with sufficient accuracy. The rapid 
determination of the bending mode frequencies is essential 
in order to place the notch filter at the proper bending 
mode frequency before the bending mode signals can degrade 
the missile attitude and rates to a point from which it 
cannot recover within the allowable flight trajectory 
envelope. The allowable time for the initial determination 
of the bending mode frequencies is a function of the initial 
frequency and bandwidth of the individual notch filters. 

In general the time available to eliminate the effects of 
the bending modes before the Trident Missile deviates from 
an acceptable flight trajectory is one-half second. The 
requirement for accuracy is essential to allow the minimum 
bandwidth notch filters to be used. As successive notch 
filters are added to the control system to attenuate the 
bending mode signals, an increasing amount of phase loss 
will be experienced at the frequencies used by the control 
system resulting in a decrease in the already minimal phase 
margin for overall missile stability. If very narrow notch 
filters are used to minimize the phase loss at the control 
system frequencies, then the center frequency of the notch 
must be placed much more accurately to be effective. 
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Four methods v>^ere investigated as to their potential to 
be used for inflight bending mode frequency identification 
and tracking: 

1. Reference [3] summarizes the use of Fourier 
integral filters to measure the power spectral density at 
integer frequencies. Second order curve fitting techniques 
were then used to estimate the bending mode frequencies 
from the output of the filters. The major drawback of this 
technique is that it required at least one-half second of 
integration time for each repeated frequency determination. 
Reference [3] did not address the ability of this technique 
to accurately track a dynamic bending mode, although the 
data indicated that a bending mode that was changing as 
little as 0.2 hertz per second would cause sufficient dis- 
tortion of the frequency determination as to render the 
technique ineffective. 

2. The author conducted a feasibility study of 
using a Fast Fourier Transform (FFT) algorithm to determine 
the spectrum of the rate gyro output signal in order to 
determine the bending mode frequencies. The method was 
disregarded for two reasons. The rate gyro signal would 
have to be sampled for ten seconds to obtain a minimum 
frequency resolution of 0.1 hertz and, the computational 
capability required to execute the FFT algorithm was felt 
to be excessive. 
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3. The feasibility of using a phaselock loop to 
lock-on and track a bending mode frequency was investigated 
by the author. Simulation of a phaselock loop with an eight- 
hertz center frequency and a three-hertz acquisition range 
revealed that while the tracking characteristics of the 
phaselock loop were very good, the ability to achieve lock 
with a random phase signal at the extremes of its acquisition 
range was unreliable if the signal was contaminated with fre- 
quency components from the plant and a second bending mode 

at anticiapted signal levels. The technique was not pursued 
further but with the addition of prefilters on the phase- 
lock loop input this technique has the potential to be more 
effective than the technique used in this thesis. Consi- 
derable refinement in the simulation of the phaselock loop 
nonlinearities would be required. 

4. The technique used in this thesis to identify 
and track the bending mode frequencies uses a microprocessor 
to compute the frequency of each bending mode by use of 
digital bandpass filters and linear estimation of the zero 
crossing times. The use of discrete signal processing by 
the microprocessor will not require additional sampling 
circuitry since the missile attitude and rate gyro signals 
are currently sampled and put in digital form in the Trident 
Missile for use in the digital autopilot. As shown in 
Figure 1-1, the microprocessor uses the sampled rate gyro 
data by passing it through a set of parallel digital 
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Figure 1-1 Microprocessor Algorithm Modules 
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bandpass filters in order to decompose the signal into fre- 
quency bands that are only wide enough to contain one 
bending mode frequency. The output of each of the bandpass 
digital filters is examined after each sample to see if the 
amplitude of the signal in that band has changed signs. If 
the amplitude has changed signs then the time of the zero 
crossing is linearly estimated based on the amplitude of the 
signal before and after the zero crossing and the sampling 
period. The half period detected is converted to a fre- 
quency and averaged with the frequency that the digital 
notch filter devoted to that band is currently tracking. 

The average of the current and computed frequency is then 
used to compute new coefficients such that the digital notch 
filter will be centered at the new average frequency. The 
details of the algorithm and design are discussed in Chapter 
Three and the concept is demonstrated by simulation in 
Chapter Four. 
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II. MODEL DEVELOPMENT 



In order to verify the ability of the tracking notch 
filter concept to stabilize the missile in the presence of 
bending mode feedback, a model of the system was developed 
to be used as the basis for simulation. 



A. THE MODEL OF THE PLANT 



Figure 2-1 is a block diagram representation of the re- 
duced order pitch plane attitude control loop of the Trident 
I Missile autopilot. The model ^^^as provided by Lockheed 
Missile and Space Company by reference [3] , and \\^as modified 
only to the extent that the signs at the summers for x^^ and 
x„ were changed to allow all the model coefficients to be 
considered positive quantities. Table 2-1 shows the time 
varying values for the parameters used in this model. The 
model shown in Figure 2-1 represents the rigid body dyna- 
mics of the missile and is referred to as the plant. When 
bending modes are added to the rigid body model the resulting 
composite model is referred to as the system. 

The state space representation of the plant for any time 



IS : 



Xi = 



^2 = 
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Figure 2-1 Pitch Plane Attitude Control Loop 
(Plant) 
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TABLE 2-1 



MISSILE PARAMETER VALUES AT VARIOUS TIMES INTO THE FLIGHT 



PARAMETER 


T=25 SEC. 


T=40 SEC. 


T=50 SEC. 


M 

a 


18.625 


28.625 


17.5 




22.51 


24.17 


32 . 0 


z 

a 


88.87 


180.0 


140.0 


7 

"'5 


119.19 


113.75 


137.5 


V 


1750 


2500 


3400 


rg 


0.06 


0.06 


0.06 


'^GA 


0.06 


0.06 


0.06 


*"8 


1.5 


1.48 


0.946 


CD 


0.2 


0.2 


0.2 



NOTE: The units of each parameter are listed on the list 

of symbols. The signs of the parameters are all 
positive to be consistent with the structure of the 
model used. 
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where : 

is the missile pitch angle. 

^2 is the missile pitch angle rate. 
Xj is the angle of attack. 



Using matrix notation the plant is then described by 



AX + B0 
^ c 



0 = CX 
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The transfer function relating the missile pitch angle 



9 to the commanded pitch angle 9^ is: 

e(s) 

= H(s) 



-K^[M.a+CM~Z -M Z,)] 

9 '■ 5 ^ 5— a a-5-^ ■' 

H(s) = - 3 —^ 2 

s +( a+K.K „M.)s^+[M +K .K-CM.Z -M Z J +K„M.] s+K„ (M ^ -M^Z .) 
^-rr- 9 r9 .5 '•a r9 9 ^ 5-a a— 5-^ 9 o-" 9^ 5— a a— 5"^ 



V 



V V 



V V 



( 1 ) 



The transfer function can be factored with a root finder 
and used to plot the poles and zeros of the transfer function 
in order to provide an intuitive feeling as to the response 
of the system not otherwise apparent with the state-space 
representation. The state-space representation of the sys- 
tem is used primarily in the simulation programs to obtain 
the time-domain solution to the differential equations of 
motion . 

B. THE CHARACTERISTICS OF THE PLANT 

The model of the plant given in Figure 2-1 is a model of 
the rigid body dynamics including two Newtonian integrators 
and the positive feedback of the angle of attack, x_. The 
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output pitch angle 6 is opposite in sign to the commanded 
pitch 9^. This reversal is compensated for in the outer 
guidance loop not included in this model. The sign of the 
forward gains and Z^, are normally specified as negative 

quantities which results in a net negative feedback. For 
computational convenience the gains are treated as positive 
quantities and the signs at the summers have been changed. 

The sign of the rate and position feedback summers remains 
positive to clearly indicate the nature of the positive 
feedback of the bending modes from the pitch rate and pitch 
angle gyros. 

'■V 

The model of the plant was designed to simulate the 
transient behavior of the system and not the steady state 
behavior. As a result the simulation studies are only valid 
for about two seconds of transient behavior. 

To establish a baseline or nominal performance the plant 
was simulated for a ramp and a step input with system para- 
meter values fixed at the time = 40 second values. The 
reasoning behind the selection of inputs and parameter values 
is discussed in Chapter Four. 

Selecting the appropriate values for the system para- 
meters from Table 2-1 for time = 40 seconds, the transfer 
function, given by equation (1) becomes: 

-35.77(5+0.01779) 

H(s) = (2) 

(S+0. 0098) (s+3. 608+/-J 7.17) 
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The poles and zeros of the transfer function are plotted 
in Figure 2-2. 

An examination of the poles and zeros of the transfer 
function indicates the system is stable with a pair of com- 
plex poles having a natural frequency = 8.02 RAD/SEC and 
a coefficient of damping ? = 0.449. There is also a very 
slow pole at s = -0.0098 and a zero at s = -0.01779. 

To establish the baseline or nominal performance of the 
plant without bending modes, the plant was simulated using 
the IBM DSL/360 language. The values of the plant para- 
meters from Table 2-1 for time = 40 seconds are used through- 
out this thesis. A detailed discussion of the simulation 
programs and the selection of the plant parameter values is 
contained in Chapter Four. Figure 2-3 illustrates the plant 
response to a step input pitch command of 0.1 radians. 

Figure 2-4 illustrates the response of the plant to a ramp 
input of 0.05 radians/sec. 

Note: It will be the convention in this thesis when 
presenting computer generated output in the form of a figure 
to include the computer program that was used to generate 
that figure in the computer program appendix properly 
labelled to indicate the figure number. If reference is 
made to the program it will be referred to as program 
(figure number) . 
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Figure 2-2 Locations of the Poles and 
Zeros of the Plant 
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System: Third order missile model without 

bending modes 

Input: 0.1 Radian step pitch command 

Figure 2-5 Missile response for a step command 
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Time, Sec. 



System: Third order missile model without 

bending modes 

Input: O.OSt radian/sec ramp pitch command 



Figure 2-4 Missile response for a ramp command 
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Although all simulations in this thesis are performed 
using the third-order model of the plant, it is possible to 
reduce the model to a second-order model without introducing 
significant error. The second-order model provides a con- 
venient preliminary design model since the characteristics 
of second-order systems are well established. The second- 
order model is derived by analysis of the time domain res- 
ponse of the third-order system and eliminating the effects 
of the slow pole and zero. The validity of the second- 
order model is established by comparison with the third-order 
model simulation results. 

To obtain the time domain response to a 0.1 radian step 
input the transfer function of equation (2) will be multi- 
plied by 0.1/s and the resulting function expanded in a 
partial fraction expansion. 



0.1 -5. 577(s+0. 01779) 

9(s) = H(s) = 

s s(s+0. 0098) (S+3.608+/-J 7.17) 



xvhich yields: 



-0.0098t -3.608t 

9(t) = -0.1+0.0453 e + 0.06218e COS(7 . 17t-0 . 453) 
at t = 2.0 sec . 

9(2) = -0.1 + 0.0453(0.98) + 0 . 06218 (0 . 00072) COS (14 . 34- 0 . 453) 



= -0.552 
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Examining the time domain response one finds that the 
influence of the slow pole is so slow that one may consider 
the response of the system due to that pole to be essentially 
a constant for the two second duration of the simulation with 
less than two percent error. Treating the term due to the 
slow pose as a constant reduces the time-domain solution to: 

- 3 . 6 08t 

9(t) = -0.0547 + 0.0621Se COS(7.17t - 0.453) (3) 

The system now can be considered as second-order and the 
response can be very closely approximated from tabulations 
of the characteristics of second order systems [4] . 

The transform of equation (3) is 

L[eft)]=6(3)= ♦ ■ (0-0 , 5 . 5 , 8 . 9) . ( . 3 . ^7 . .098) , 

^ s" + 7.226s * 64.4265 

Placing all terms over a common denominator: 

^ -0.0547 s“-0.3952 s-3.524 + 0.05589 s“ +0.3967 s 

0(S) = 

s (s“+7.226 s+64.4265) 

2 

Cancelling the essentially equal s“ and s terms in the 
numerator and factoring out the transform of the 0.1 radian 
step input (0.1/s) yields the following reduced order trans- 
fer function relating the missile pitch angle 9 to the 

commanded pitch angle 9 . 

o c 



= H (s) = 

^ s“ + 7.226 s + 64.4265 



(3a) 
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The steady state gain of the system is: 



s H^(s) 



LIM 

§■^0 



=-0.547 



The peak overshoot \\rill be: 

TTC 

St = 1 * ® /T~? = 

The time of peak overshoot xvfill be 



t = 
P 



TT 



W J \ - t} 
n 



= 0.438 sec 



Using the approximation for the settling time 



. 4.0 



tW 

n 



= 1.11 sec 



The validity of the second-order model for a ramp input 
is established by comparison of the value of the missile 
pitch angle determined analytically from the second-order 
transfer function at time = 2 seconds with the value resulting 
from the simulation of the third-order plant. 

Given a transfer function in the form: 



G(s) = 



B. + B,s + B^s^ + + B^s"' 

0 X ^ in 



Aq + AjS + AjS' 



+ A s 
n 



n 
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The steady state time domain response o£ the output to a 
ramp input, u(t) = at, can be determined from the following 
formula [5] : 



y.Ct) = 



B, 



A, 



at + 



(A„Bj - 



a; 



Using the second-order transfer function determined in 
equation (3a) one finds that the output at time=2 for a ramp 
input, u(t) = O.OSt is: 

-3S.24 (0.05)t 0. 05(35. 24) (7. 226) 

e (t) = + 5 = -0.05242 

^ 64.4265 (64.4265) 



Table 2-2 summarizes the results obtained from the second 
order model with those resulting from simulation of the third 
order model. Again, it must be noted that the simulation 
studies contained in this thesis were performed using the 
third-order model of the plant. The second-order model was 
derived only as a convenient preliminary design tool for any 
subsequent studies. 
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TABLE 2-2 



COMPARISON OF THIRD-ORDER AND SECOND-ORDER MODEL RESPONSES 



CHARACTERISTIC 


SECOND- 

ORDER 

MODEL 


THIRD- 

ORDER 

SIMULATION 


MAXIMUM 

OVERSHOOT 


1.21 


1.21 


TIME AT MAX. 
OVERSHOOT 


0.44 SEC 


0.44 SEC 


STEADY STATE 
GAIN 


-0.5625 


-0.547 


SETTLING 
TIME (2%) 


1.11 SEC 


1.11 SEC 


RAMP OUTPUT 
AT T=2 SEC 


-0.0524 


-0.0526 
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C. THE MODEL OF THE BENDING MODES 



Reference [3] provided the model of the plant with one 
bending mode shown in Figure 2-5. The bending mode attitude 
signal and bending mode rate signal are added to the 
rigid body attitude signal and rigid body rate signal X 2 
respectively. The composite attitude and rate signals are 
each fed back to the system input in a positive sense and 
the overall system becomes highly unstable. 

Reference [2] states that the bending mode contamination 
of the rigid body rate signal has a much greater destabilizing 
effect than does the bending mode contamination of the rigid 
body attitude signal. The bending mode model is therefore 
based on modeling the bending mode rate and the bending mode 
attitude signal is formed by integrating the rate signal. 

Reference [6] contains the bending mode frequency res- 
ponse for a Poseidon Missile when the bending modes were 
excited by slewing the first stage TVC nozzles. Although 
the data was measured on the Poseidon Missile ;vith various 
flight configurations, the results are consistent with the 
general distribution of bending modes shown in Figure 5.5-1 
of Ref. [2]. The distribution of the bending modes indicates 
that the bandwidth of the bending modes increases as the 
center frequency of the bending mode increases. The distri- 
bution of the bending modes used in this thesis is shown in 
Figure 2-6. The low frequency bending mode may be centered 
anywhere in a band from 5 to 8 hertz and the high frequency 
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Kro 



Figure 2-5 Model of the Plant with 
One Bending Mode 
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Frequency Hz 
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Figure 2^-6 Dynamic Bending Mode Distribution 
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bending mode may be centered anywhere in the band from 8 to 
15 hertz. The high frequency bending mode has approximately 
twice the bandwidth of the low frequency bending mode. The 
simulation studies contained in Chapter Four demonstrate the 
ability of the adaptive notch filtering scheme to identify, 
track and filter bending modes that have time varying gains, 
bandwidths and center frequencies. 

The bending mode rates were modeled as the outputs of 
bandpass filters with a very high Q. The transfer function 
relating the output Y(s) to the input U(s) for a general 
second-order bandpass filter is: 

Bs 

G(s) = 

2 2 
s + Bs + W 

c 



where : 

= the center frequency of the filter 
B = the 5 dB bandwidth of the filter 

The block diagram of G(s) is shown in Figure 2-7. Y(s) 
represents the bending mode rate output and since the bending 
mode attitude output is the integral of the bending mode rate, 
an output of the bending mode attitude is also available from 
the rate model as indicated in Figure 2-8 where Y^(s) is the 
bending mode rate output and Y 2 (s) is the bending mode atti- 
tude output. 



36 




Figure 2-7 Bending Mode Rate Model 



Y.(S) 




Figure 2-8 Bending Mode Rate and Attitude Model 
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Inspection of Figure 2-8 reveals that modeling the 
bending mode rate as a bandpass filter is equivalent to the 
bending mode model shown in Figure 2-5. 

Figure 2-9 shows the system model consisting of the plant 
and t\vo bending modes that are used for the simulation 
studies in this thesis. 

D. CHARACTERISTICS OF THE PLANT WITH TWO BENDING MODES 

Preliminary values for the bending mode parameters of gain, 
center frequency and bandwidth were selected for the purposes 
of analysis and to develop some insight as to the elements 
controlling the behavior of the system. The bending mode 
parameters are varied during the simulation studies presen- 
ted in Chapter Four but the detailed analysis is not repeated 
for each parameter study. 

The bending mode parameters selected for analysis of the 
system are; 

= 20.0 

=1.10 RAD/SEC 

Wj^ = 62.83 R.^D/SEC (10 Hz center frequency) 

= 240.0 

B^ = 0.1 RAD/SEC 

W^ = 37.7 RAD/SEC (6 Hz center frequency) 

The frequency response of the transfer function relating 
the missile pitch angle 9(s) to the input pitch command 9^(s) 
with the rate and attitude feedback loops open was determined 
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Figure 2-9 Model of Plant with Two Bending Modes 



40 




for various configurations of the plant and bending modes 
shown in Figure 2-9. The plant parameters from Table 2-1 
for time = 40 seconds \tfere used. 

Figure 2-10 shows the open loop frequency response of the 
plant with no bending modes. 

Figure 2-11 shows the open loop frequency response of the 

6.0 Hz bending mode alone. 

Figure 2-12 shows the open loop frequency response of the 
plant and the 6.0 Hz bending mode. 

Figure 2-13 shows the open loop frequency response of the 

10.0 Hz bending mode alone. 

Figure 2-14 shows the open loop frequency response of the 
plant and the 10.0 Hz bending mode. 

Figure 2-15 shows the open loop frequency response of the 
plant with both the 6.0 Hz and the 10.0 Hz bending modes. 

The matri:x representation of the state equations derived 
from the closed loop system shown in Figure 2-9 is: 

« * 

X = Ax + Be 

c 

9 = Cx 
\\^here : 

is the rigid body pitch angle of the missile 
x^ is the rigid body pitch rate of the missile 
Xj is the missile angle of attack 

X is the pitch angle resulting from the lower frequency 
bending mode 
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Frequency, Hz (linear scale) 



Figure 2^-10 Open Loop Frequency Response 
of the Plant 
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Figure 2-11 



Open Loop Frequency Response Of 
Low Frequency Bending Mode 
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gure 2-12 Open Loop Frequency Response of 
the Plant and Low Frequency 
Bending Mode 
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Figure 2-13 Open Loop Frequency Response of 
the High Frequency Bending Mode 
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dB 




Figure 2^14 Open Loop Frequency Response o£ 
the Plant and High Frequency 
Bending Mode 
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dB 




Figure 2^15 



Open Loop Frequency Response of 
Plant with two Bending Modes 
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X|- is the pitch rate resulting from the lower frequency 
^ bending mode 

X, is the pitch angle resulting from the higher fre- 
quency bending mode 

x.^ is the pitch rate resulting from the higher frequency 
bending mode 

9^ is the commanded pitch angle 
0 is the missile pitch angle 

The A matrix is ; 



All 


A12 


A13 


A14 


A15 


A16 


A17 


A21 


A2 2 


A23 


A24 


A2 5 


A26 


A27 


A31 


A3 2 


A33 


A3 4 


A35 


A3 6 


A37 


A41 


A42 


A43 


A44 


A4 5 


A46 


A4 7 


A51 


A52 


A53 


A54 


ASS 


AS 6 


AS7 


A61 


A62 


A63 


A64 


A6S 


A66 


A67 


A71 


A72 


A73 


A74 


A7S 


A7 6 


A77 



The B matrix is : 

[B1 B2 B3 B4 B5 B6 B7]^ 

The C matrix is : 

[Cl C2 C3 C4 C5 C6 C7] 
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where : 



All = 

A12 = 

A13 = 
A14 = 
A15 = 
A16 = 
A17 = 

A41 = 

A42 = 

A43 = 
A44 = 

A45 = 

A46 = 
A47 = 

A71 = 
A72 = 
A73 = 



0.0 

1.0 

0.0 

0.0 

0.0 

0.0 

0.0 



A21 

A22 

A23 

A24 

A25 

A26 

A27 









= M 



a 



■’^GA^e^^d 



-Ip 

^rG r9 6 5 

-'^GA‘6‘''5 

-Z-M. 
^rG r6 6 o. 



A31 = 

A3 2 = 

A33 = 
A34 = 
A35 = 
A3 6 = 
A3 7 = 



K 

®TT 



z 



V 

■'^GA^eir 

I 7 -7 

■'^rG'^r9"0 V 

f - ^ 

■’^GA“'9 V 

n 

^ 5 

■’^rG^r9^elT 



0.0 


A51 


°L®L 

’^GA 

^9 

= GtBtK 


A61 =0.0 


0.0 


A52 


A62 = 0.0 


0.0 


A53 


L L 
= 0.0 


A64 = 0.0 


0.0 


A54 


= 


A6S = 0.0 


1.0 


ASS 


= G, B,ij; 


-B, A66 = 


0.0 


AS6 


L L^rG 
= G, B.K^ 


L 

A66 = 0.0 


0.0 


AS7 


L L 9 ^ 

= G, B, K „K. , 


A67 = 1.0 




L L r9 


K- 


°»®h*ga 


B1 = 


O 

O 


Cl = 1.0 


K. 




B2 = 




C2 = 0.0 


0.0 


B3 = 


Z. 

- 


C3 = 0.0 
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A 74 = 


'^-ra 


B 4 = 0.0 

^0 

B 5 = G.B,^ 

L L'^ga 


C 4 = 

rg 


A 75 = 




C 5 = 0.0 


A 76 = 


Vnh-%' 


B 6 = 0.0 


= ^GA 


A 7 7 = 




^9 

B 7 = G B ^ 
” ”*GA 


C 7 = 0.0 



The closed loop transfer function H(s) relating the 
missile pitch angle 9 to the commanded pitch angle 9^ is: 



H(S) = 



N’ s^ + N.s^ + N-s-^ + N.s^ + N^s + N. 

5 4 ^ 2 1 0 

s^ + D.s^ +D.s^ + D.s^ + D^s'^ + D„'s^ + D, s + D„ 
DO 4 o 2 10 



where : 

? ^ ^ 7 ? 

N- = -M„W„^-M,B, B„+B„(M ^-M.-#)+B-(M ~M.~) -M .W. +G, W “B, 
3 5H 5LH H^aV 5V L aV SV 5L LH L 



N. = 



Z Z ^ 



2,,. ^3 
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7 

2.,. ^5: 



+W, G„B„M 
L H H a 



N = W ^ CM — -M — ) 
‘^0 H L '• a V S'V '* 



D = B, + B„ + ^ - K^K ^N 



H V 



D. = W, "' + B,-t^+M +B, B, +B„~+W,,‘'-K„K ^N,-K^N 



re* 5 

z 

a 



D. = 



'L LV a L H HV H 9 r0 4 9* 5 






2 a.,., 2,., 2 



W-^-M +B„B,M +BuW, "+W„“B,-^+W„“M -K,K .N.,-K.N, 
LaHLa HLV H}j[ H LV H a 9r9 2 9a 



D- = B„W, +B,W„^M +W, ^W„^^-K,K -K^N- 

2 KLaLHaLHV 9r9 192 



= VSk-''e''reNo-KeNi 



Do = -Vo 



For the bending mode parameters selected and the plant 
parameters for time = 40 seconds the transfer function H(s) 
is ; 



_34.97 (s+4 8.4) (s- 48.7) (s+ 0.018) (s-^0. 215 + /- j 48. 65) 

'' ^ 9^ (Sj (s-^O. 0iy(s+3.63 + /- j 6.997) (s-3 . 54-^/- j38 . 15) 

^ (S-2.91+/-J62.25) 
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The transfer function H(s) Vvas factored using the "roots” 
program contained in Ref. [7] . The location of the poles was 
confirmed by finding the eigenvalues of the A matri:x using 
the "BASMAT" program also contained in Ref. [7] . The poles 
and zeros of H(s) are plotted in Figure 2-16. 

The system composed of the plant plus bending modes is 
clearly unstable with four poles of the closed loop transfer 
function in the right half of the S-plane. 

Figure 2-17 shows the simulation the system shown in 
Figure 2-9 for a step input of 0.1 radians. Line 1 of 
Figure 2-17 is the response of the plant without bending 
modes and line 2 is the response of the system composed of 
the plant with bending modes. The vertical scaling in Figure 
2-17 was restricted to keep the response of the plant in 
proper perspective. The restriction of the vertical scale 
causes saturation of the plot of the system response (line 2). 
Figure 2-18 shows the response of the system with the verti- 
cal scale extended to remove the effects of saturation of 
the plot. For clarity, the response of the plant without 
bending modes was not included in Figure 2-18. 

To validate the simulation, the time response of the 
state equations with a 0.1 radian step input was determined 
by using the "RTRESP" program contained in Ref. [7] . The 
pitch angle of the missile, 0 was found to be described by 
the equation; 
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Figure 2^16 Location of Poise and Zeros 
of the Plant with Two 
Bending Modes 
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System; Third order missile model with two 
bending modes and no adaptive notch 
filters 

Input; 0.1 radian step command 

Line 1 is the nominal missile response without 
bending modes 

Line 2 is the missile response v\:ith two bending 
modes 

Figure 2rl7 Missile Response with Two Bending Modes 
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System; Third order missile model with two 
bending modes and no adaptive notch 
filters 

Input: 0.1 radian step command 



Figure 2-18 Missile Response with Two Bending 
Modes (Expanded Scale) 
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0(t3 = -0.09952 + 0.0457 e * 

+ 0.0582 cos (6.997t-0.41) 

+ 0.00246 cos (38 . 15t- 2 . 45) 

+ 0.00093 cos(62.25t+3.06) 

At time = 1.0 seconds, the above equation indicates the 
pitch angle will be 9(1) = -0.10162 radians. The numerical 
results of the simulation used to generate Figure 2-17 indi- 
cate the pitch angle to be 9(1) = -0.1017 radians. The 
proceeding results indicate the simulation program is valid. 

The system was simulated with only the high frequency 
bending mode and with only the low frequency bending mode. 

The purpose of these simulations was to shov\: that each bending 
mode would by itself destabilize the system and that the fre- 
quency of the oscillation of the resulting unstable system 
would reflect very nearly the center frequency of the 
bending mode that destabilized the system. Figure 2-19 illus- 
trates the response of the system for a 0.1 radian step 
input when only the high frequency bending mode was present. 
Figure 2-20 shows the same simulation with an e:xpanded verti- 
cal scale. Figure 2-21 shows the response of the system for 
a 0.1 radian step input with only the low frequency bending 
mode present. Figure 2-22 shows the same simulation with 
an expanded vertical scale. Figure 2-23 shows the system 
response \vith both bending modes for a ramp input 
9(.Ct) = 0.05t. 
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Time, Sec. 

System: Third order missile model with 10 Hz 

bending mode and no adaptive notch 
filter 

Input: 0.1 radian step command 

Line 1 is the nominal missile response without 
bending modes 

Line 2 is the missile response with 10 Hz bending 
modes 



Figure 2-19 Missile Response with 10 Hz Bending Modes 
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filter 

Input: 0,1 radian step command 



Figure 2-20 Missile Response with 10 Hz Bending 
Mode (expanded scale) 
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System: Third order missile model with 6 Hz 

bending mode and no adaptive notch 
filters 

Input: 0,1 radian step command 

Line 1 is the nominal missile response without 
bending modes 

Line 2 is the missile response with 6 Hz bending 
modes 

Figure 2-21 Missile Response \vith 6 Hz Bending Mode 
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System: Third order missile model x\rith 6 Hz 

bending mode and no adaptive notch 
filters 

Input: 0.1 radian step command 



Figure 2r22 Missile Response with 6 Hz Bending 
Mode (expanded scale) 
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System: Third order missile model with 6 and 10 

Hz bending modes and no adaptive notch 
filters 

Input: 0,05 radian/sec ramp command 

Line 1 is the nominal missile response without 
bending modes 

Line 2 is the missile response with bending modes 

Figure 2-23 Missile Response for a Ramp Input \vith 
Bending Modes 
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Figures 2-19 through 2-22 clearly indicate that each 
bending mode is in itself destabilizing. Analysis of the 
numerical data associated \vith Figures 2-19 through 2-22 
confirms that the frequency of the oscillations is essen- 
tially at the frequency of the bending mode center frequency. 
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III. ADAPTIVE DIGITAL FILTER DESIGN 



A. FACTORS AFFECTING THE DESIGN 

The emphasis of the design was simplicity. Significant 
improvement in the results could have been achieved by use 
of higher order filters and by using more sophisticated fil- 
ters, but the improved performance would have been accom- 
plished at the expense of obscuring the very simplicity ;vhich 
makes the proposed technique so appealing. 

The sampling frequency used throughout this thesis is 100 
hertz. This selection is dictated by the sampling frequency 
used by the autopilot in the Trident Missile. Higher sam- 
pling frequencies were found not to significantly improve 
the performance of the algorithm, however, sampling fre- 
quencies lower than 60 Hz degraded the algorithm to where 
it would not perform reliably. 

The bending mode frequencies \^ere selected such that one 
and only one bending mode would be in each of the bands of 
5-8 hertz and 8-12 hertz as shown in Figure 2-6. A priori 
knowledge of the band or range of frequencies over which a 
given bending mode may vary during the flight of the missile 
is essential in this design. If specific events during the 
flight such as staging change the bands of the bending modes 
then the algorithm may also be modified to adapt to the new 
distribution with the only requirer.ient being that the times 
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or events which change the band limitations be preprogrammed 
in the microprocessor algorithm. 

B. DESIGN OF THE ALGORITHM MODULES 

Figure 3-1 shows the' modular nature of the algorithm for 
the determination, tracking and notch filtering of the bending 
mode signals. The sequence of modules associated with the 
determination, tracking and notch filtering for each of the 
bending modes are collectively termed a channel. The channels 
are identical except for the numerical values of the coeffi- 
cients; therefore, only the design of the 8-12 hertz (high) 
channel is presented in detail. The results of the design 
of each of the modules comprising the low channel are pre- 
sented when different from the corresponding module in the 
high channel. 

1 . The Bandpass Filter Module 

The overall bandpass filter was designed as two 
cascaded second-order bandpass sections. The first section 
is a first-order Butterworth lowpass Filter transformed into 
a bandpass filter and then transformed into the z-domain 
using the bilinear transformation. The second section is 
also a bandpass filter designed to adjust the overall shape 
of the gain curve to the desired form. 

A first-order Butterworth loivpass Filter normalized 
to a cutoff frequency of 1 radian/second has a transfer 
function ; 
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Figure 3-1 Adaptive Notch Filter Algorithm Modules 
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HipCs) 



1 

s + 1 



The transformation from lowpass to bandpass is: 



s 






2 

s 



+ W ^ 
c 

Bs 



where : 

B is the 3 dB bandwidth defined as W - W. 

c L 

W is the center frequency defined as / W W. 

is the upper 5 dB frequency in rad/sec and is the 
lower 3 dB frequency in rad/sec. 

The resulting analog domain transfer function for 
the bandpass filter is: 



HbpCs) 



— 
s + 



Bs 

7 

Bs + W “ 
c 



(4) 



Before applying the bilinear transformation to trans- 
form the analog bandpass filter to the discrete domain it 
is necessary to prewarp the center frequency [8] . The 
prewarping is done by using the following transformation; 



2 TAN 



“a = 



(^) 



(5) 
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where : 



is the desired discrete or digital frequency defined 

as : 

T is the sampling frequency. 

1/2 

W is the analog center frequency defined as W = (W ) 

Selecting 8.0 hertz as the lower cutoff frequency and 
12.0 hertz as the upper cutoff frequency. The values to be 
used in equations (4) and (5) become: 

Wj^ = 2tt* 8.0 = 50.265 rad/sec 

W = 2tt* 12.0 = 75.398 rad/sec 
u 

W = (W W, = 61.652 rad/sec 
c u L 

B = W -W^ = 25.133 rad/sec 
u L 

Applying equation (5) with the selected value of 
yield the prewarped center frequency: 

= 63.582 rad/sec 

Substituting the above value, equation (4) becomes: 

H(s) = ^ (6) 

s + 25.133 s + 4042.77 
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The bilinear transformation for the s-domain to the 



z-domain is; 



- 2(Z - 1) 

s T(Z + 1) 



Applying the transformation to equation (6) , the dis 
Crete transfer function for the bandpass filter becomes: 

H (Z) = 0.105(Z^ - 1) 

^ Z^-1.4677 Z + 0.79405 

After slight adjustments to obtain the desired gain 
and crossover frequency with the low frequency channel band- 
pass filter (see Figure 3-9), the transfer function H^(Z) is 






0.1 (Z - 1) 



1.45625 Z + 0.81 



The pole and zero configuration of H^(Z) is shown in 
Figure 3-2 and the frequency response of H^(Z) is shown in 
Figure 3-3. 

As expected, the frequency response of H^(Z) is too 
flat to provide the needed stopband attenuation for the 
adjacent channels. The frequency response of the overall 
bandpass filter is greatly improved by adding a second sec- 
tion of filter with a sharp resonant peak at the center fre- 
quency. A bandpass filter with the desired resonant peak 
was designed directly in the z-domain using a digital filter 
frequency response program (program 3-5) to adjust the poles 
and zeros to obtain the desired characteristics. The trans- 
fer function of the second filter section is: 

H (Z) = 0-456(Z^ - 0.809 Z + 0.25) 

^ Z^ - 1.3753 Z + 0.7225 
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Figure 3-2 Discrete Bandpass Filter 
Pole/Zero Locations 
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Figure 3-3 Frequency Response of the First Section 
of the High Frequency Digital Bandpass 
Filter 
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The pole and zero configuration of H 2 CZ) is shown in 
Figure 3-4. The frequency response of H^C^) is shown in 
Figure 3-5. 

Figure 3-6 shows the flow graph of the algorithm for 
the implementation of the two cascaded sections and Figure 
3-7 shows the resulting frequency response of the cascaded 
sections . 

Figure 3-8 shows the frequency response of the low 
channel bandpass filter designed in the same manner. Figure 
3-9 illustrates the resulting band coverage of each channel. 

The possibility of a bending mode occurring at the 
eight-hertz cross-over frequency between channels causing 
both channels to lock on that frequency and leaving one bend- 
ing mode unattenuated is compensated for by initializing both 
notch filters at the cross-over frequency, thereby attenuat- 
ing a bending mode at eight-hertz immediately and freeing 
either channel to track and remove the unattenuated bending 
mode . 

2 . The Zero Crossing Detector Module 

The zero crossing detector is based on a linear 
estimate of the crossing time if the amplitude of the output 
of the bandpass filter changes signs between samples. Linear 
estimation of the crossing time is a feature incorporated 
into the digital simulation language program (DSL) as the 
cross function but could also be programmed using the 
algorithm shown in Figure 3-10. The potential for using 
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UNIT CIRClU 




Figure 3^.4 Discrete Bandpass Filter 
Pole/Zero Locations 
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Figure 3-5 Frequency Response of the Second Section of 
the High Frequency Digital Bandpass Filter 
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Figure 3^6 Digital Filter Implementation 
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Figure 3-7 Frequency Response of Cascaded Filter 

Sections of the High Frequency Bandpass 
Filter 
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Figure 3-8 Frequency Response o£ Low Frequency 
Digital Bandpass Filter 



76 




§ 



0.00 sToO I'.OO s'. 00 s'. 00 lb. 00 12.00 ik.oo To. 00 

Frequency, Hz (linear scale) 



Figure 3-9 Frequency Coverage of High and Low 
Frequency Digital Bandpass Filter 
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Figure 3-10 Zero Crossing Time Algorithm 
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a sinusoidal curve estimation rather than a linear estima- 
tion was considered and dismissed because any increase in 
the accuracy of the estimate of the zero crossing time was 
offset by the jitter in the actual crossing time caused 
by leakage into the channel of the adjacent bending mode and 
plant frequencies. 

The major deficiency in the use of a zero crossing 
detector as a frequency measuring device is the nonlinear 
error in the frequency estimate as a function of the actual 
frequency. For example, if the frequency being measured is 
5.0 hertz the half-period measured should be 0.1 seconds. 

If noise caused a 0.01 second measurement error and the 
measured half-period was 0.09 seconds then the frequency 
detector would estimate the input frequency to be 5.55 hertz. 
If the input frequency was 10.0 hertz the same 0.01 second 
error would result in an estimate of the input frequency of 
12.5 hertz, a five times greater error for the same measure- 
ment error. 

The effect of the nonlinear quality of the estimated 
frequency in the presence of noise causes a much greater 
variation in the estimate of the higher frequency bending 
modes and contributes to the limit cycle problem described 
below. 

3 . The Frequency Averaging Module 

The primary source of limit cycle behavior in the 
tracking algorithm appears to be due to leakage of an 
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adjacent bending mode frequency into the zero crossing detec- 
tor causing a periodic error in the frequency determination. 
The frequency averaging module simply determines the mean 
of the last detected frequency from the frequency detector 
and the current notch filter center frequency. The mean of 
the two frequencies becomes the new notch filter center fre- 
quency. The process of averaging the last detected frequency 
with the current notch filter center frequency reduces the 
amplitude of the limit cycle but does not eliminate it. To 
further reduce the amplitude of the limit cycle, several 
approaches may be taken; 

(a) Greatly improve the stopband attenuation and roll- 
off characteristics of the bandpass filters. This approach 
is discussed in the conclusions contained in Chapter Five. 

(b) Average the detected frequency over a sufficient 
number of zero crossings before moving the notch. This 
approach was found to significantly reduce the amplitude of 
the limit cycle but was so slow in moving the notch to the 
correct frequency that the bending modes would increase to 
such an amplitude that the missile \>rould be destroyed 

4 . Design of the Notch Filter Module 

The design of a discrete notch filter and the calcu- 
lation of the coefficients is simple and efficient. Figure 
3-11 illustrates the pole and zero location of a second- 
order notch filter on the Z-plane. Since the zero is 
located on the unit circle, the center frequency of the 
notch is completely specified by the angle 
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FILTEI1 POLE ZERO LOCfiTIONS 



Figure 3-11 General Pole/Zero Locations 
for a Second-Order Discrete 
Notch Filter 
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The angle is referred to as the discrete fre- 
quency and is defined as the angle from the positive real 
axis to the zero such that: 



Q = W T 
c 



where : 

is the desired notch center frequency in rad/sec 
T is the sampling period in seconds 

The location of the zero on the unit circle at a 
selected discrete frequency, 9.^, ensures that the center 
of the notch will be at the desired frequency completely 
independent of the selected location of the poles. The loca 
tion of the poles on the same radial (discrete frequency) 
as the zero maintains the approximate s>Tnmetry of the notch. 
The radius of the poles determines the bandwidth of the 
notch. The closer the poles are to the unit circle, the 
narrower the notch will be and the longer the settling time. 

The transfer function H,,(Z) for the notch filter is 

N 

of the form: 



Hj^(Z) 




M Z + 1) 
B1 Z + B2 



where the coefficients are defined by the following formulas 

A1 = 2 cos (W^T) 

Bl = 2 r cos (W T) 

p ^ c ^ 

B2 = r - 

P 
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r is the radius of the poles 
P 

is the desired center frequency of the notch in rad/sec 
G is the gain 

For unity DC gain H(Z)jZ=+l=l yields the formula 

. 1 - B1 + B2 

^ 2 - A1 

The relationship between the pole radius and the 3 dB 
bandwidth of the notch filter is a function of center fre- 
quency but can be approximated by making certain simplifying 
assumptions. Figure 3-12 shows the geometry on the unit 
circle when the discrete frequency, is at the 3 dB 

cutoff frequency. The gain at the cutoff frequency is: 



G = 




0.707 



( 7 ) 



In this application it is reasonable to assume that 
jZ^l^jp^l, the poles are close to the unit circle, and the 
bandwidth of the filter is small (W^TlW^T) . With the above 
assumptions, equation (7) then reduces to: 



G = 




0.707 



C8) 



2 

Using the pythagorean theorem and substituting Z^ = 

2 

0.5 from equation (8) the following relationship can be 

found : 
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Figure 3-12 Vectors Used to Calculate 
Notch Filter Bandwidth 
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Cl - rp) = Z1 



C9J 



For small angles the cord can be assumed to be equal 
to the ARC length, therefore: 

Z1 = W^T - WjjT = 0.5 Bjj = 0.5 (10) 

where : 

Bj^ is the discrete bandwidth in rad 
B^^ is the analog bandwidth in rad/sec 

Substituting equation (10) into equation (9) and 
solving for B^ yields the approximate relationship between 
the notch filter analog frequency bandwidth and the pole 
radius in the discrete domain. 

2(1 - r ) 

®W ° T-^ 

The numerical data associated with the frequency 
response curves for the notch filters indicates that equation 
(11) yields bandwidths that are larger than the desired 
bandwidths but the error is less than 10 percent in the 
frequency range of 5 to 15 hertz. 

Figure 3-13 shows the gain vs frequency response of 
the low channel notch filter with a center frequency of 6.0 
hertz and a pole radius of 0.9. This configuration has a 
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Figure 3-13 Frequency Response of Digital Notch 
Filter with 0,9 Pole Radius Centered 
at 6 Hz 
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3 dB bandwidth of 2.95 hertz. Figure 3-14 sho^^^s the phase 
vs frequency response of the same notch filter. 

Figure 3-15 shows the gain vs frequency response of 
the high channel notch filter \vith a center frequency of 
10.0 hertz and a pole radius of 0.8. This configuration has 
a 3 dB bandwidth of 5.9 hertz. Figure 3-16 shows the phase 
vs frequency response of the same notch. 

The argument used to create Figure 2-6 established 
that the higher of the tivo bending modes \>fould have a vifider 
bandwidth than the lower frequency bending mode, therefore 
the bandwidth of the notch filter for the high channel was 
designed to be wider than the low channel. There is a limit 
to which the notch bandwidth can be increased without in- 
ducing excessive phase loss at the plant operating frequency 
of two-hertz. For the notch bandwidths used in this algorithm 
the high channel, with the notch centered at 10.0 hertz, in- 
troduces a phase loss of 7.75 degrees at two-hertz and the 
low channel, with the notch centered at 6.0 hertz, intro- 
duces a phase loss of 11.0 degrees at two-hertz. The com- 
bined phase loss introduced by the notch filters used in 
this algorithm does not have a pronounced effect on the 
stability of the plant. 

5 . Flow Chart of the Algorithm 

Figure 3-17 provides a flow chart of the entire 
algorithm. The "prime time" computations are the processing 
of the error signal and outputting the filtered value to the 



87 




Frequency, Hz (linear scale) 



Figure 3-'14 Frequency Response of Digital Notch 
Filter with 0,9 Pole Radius Centered 
at 6 Hz 
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Figure 3-15 Frequency Response of Digital Notch 
Filter with 0.8 Pole Radius Centered 
at 10 Hz 
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Frequency, Hz (linear scale) 



Figure 3-16 Frequency Response of Digital Notch 
Filter with 0.8 Pole Radius Centered 
at 10 Hz 
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Figure 3-17 Algorithm Flowchart 
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autopilot with minimum processing delay. The "slack time" 
computations are the computations performed after the error 
signal has been processed and which do not introduce pro- 
cessing delay in the control system. With a 100 hertz 
sampling frequency, sufficient time should be available to 
complete all slack time calculations every cycle, with 
currently available microprocessors. 

C. TEST OF THE ALGORITHM 

The algorithm was first tested as shown in Figure 3-18. 
An input signal consisting of the sum of eqnal amplitude 
2.0, 6.0 and 10.0 hertz sine waves initially in phase was 
used to simulate the rate gyro output signal with bending 
mode contamination. The signal was first passed through 
the tracking notch filter modules and the result was then 
used as an input to the two channels of the bending mode 
frequency determination, tracking and filtering algorithm. 
The unusual process of passing the input signal through the 
notch filters first ivas done to allow the effects of feed- 
back to be observed as the notch filters remove the very 
signals that the remainder of the algorithm is trying to 
track. 

Figure 3-19 shows the two-hertz input signal (line 3) 
to the notch filters and the resulting output signal (line 
2, accentuated). The notch filters \\fere reasonably able 
to filter out the 6 and 10 hertz components and pass the 



92 



I 




10 Hz 



Figure 3-18 Adaptive Notch Filter 

Algorithm Test Arrangement 
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Line 1 is zero amplitude reference 

Line 2 is the output of the adaptive digital 
notch filters 

Line 3 is 2 Hz input signal 



Figure 3-19 Output of the Adaptive Digital Notch 
Filters 
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two-hertz. The numerical data associated with Figure 3-19 
indicates an average phase loss of 18 degrees was experienced 
over the 2 seconds of simulation. 

Figure 3-20 shows the tracking history of each of the 
channels. The low channel locked on the six-hertz frequency 
(actual lock on frequency was 6.0546 Hz) after 0.56 seconds 
and remained locked on for the remainder of the simulation. 

The high channel went into a four cycle per second limit 
cycle with an average frequency of 9.85 hertz. The cause 
of the limit cycle is leadage of the two-hertz signal when 
the amplitude of the two-hertz reaches near maximum ampli- 
tudes. The limit cycle has a minimum effect on the output 
of the notch filters since the high channel notch bandi^idth 
is nearly six-hertz and the magnitude of the limit cycle 
is only 1.65 hertz. 

Figure 3-21 shows the input (line 1) and output (line 2, 
accentuated) of the low frequency channel bandpass filter. 
Figure 3-22 shoir^s the input (line 1) and the output (line 2, 
accentuated) of the high channel bandpass filter. The effects 
of the leadage of the two-hertz signal can be seen in the 
output of both filters. 
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Figure 3^20 Frequency Detector Output for 6 and 
10 Hz Input 
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Line 1 is the output from the adaptive notch 
filter 

Line 2 is the output from the low channel 
bandpass filter 



Figure 3-21 Low Channel Bandpass Filter Output 
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Line 1 is the output from the adaptive notch 
filters 



Line 2 is the output from the high channel 
bandpass filter 



Figure 3-22 Output from the High Channel Bandpass 
Filter 
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IV. SIMULATION STUDIES 



A. DIGITAL SIMULATION LANGUAGE (DSL/360) 

The programming languages used throughout this thesis are 
the IBM Digital Simulation Language (DSL/360) [9] and FORTRAN 
IV. DSL/360 has not been released for general use by the IBM 
Corporation but was used in this thesis for the simulation 
studies due to its enhanced capabilities. The exclusive 
features of DSL/360 were used only in modeling of the con- 
tinuous systems. The modules of the adaptive digital filter 
are programmed exclusively in FORTRAN IV and may be used with 
any continuous systems simulation language that will accept 
FORTRAN subroutines. 

B. BLOCK DIAGRAM OF THE SIMULATION MODEL 

Figure 4-1 illustrates the block diagram of the system 
used to develop Program 4-2. The block diagram forms the 
core program for all the simulation studies in this chapter 
with the only changes being the addition of a ramp function 
for the input command in Program 4-3 and ramp functions to 
create dynamic bending mode frequencies in Programs 4-10 and 
4-12. 

C. SIMULATION CONDITIONS 

The parameters of the plant selected for all the following 
simulation studies are from Table 2-1 at time = 40 seconds. 
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Figure 4^1 Simulation Study Model 
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The selected plant parameters represent the values when the 
missile passes through the point of maximum dynamic pres- 
sure (max Q) and is the least stable even v\rithout the pre- 
sence of bending modes. When simulation studies were 
conducted with plant parameters yielding a more stable plant 
the results indicated the proposed adaptive digital filter- 
ing technique functioned essentially the same and therefore 
only the worst case conditions are presented. 

Tvvo pitch commands from the guidance system were selected 
for simulation: a 0.1 radian step input and a 0.05 radian 

per second ramp. The step command simulates a command from 
the autopilot which would normally occur at water exit, 
staging, or inflight restart of the guidance system after a 
failure. The simulation studies with a step input are in- 
tended to simulate the condition of a guidance system 
restart at the point in the flight when the missile has the 
least stability, a highly unlikely but undoubtably severe 
condition. The ramp input simulates a more normal condition 
where the autopilot responds to \>iind shear at the point of 
least stability. The ramp input response serves to demon- 
strate that although the bending modes are not excited as 
violently as they are with a step input, the adaptive digital 
filters successfully detect their presence and provide the 
proper filtering. The ramp input again includes the simu- 
lation of an inflight restart of the guidance system since 
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the simulation is started with no history available to the 
adaptive digital filters or frequency detector algorithms. 

The bending mode frequencies were the primary variable 
in the simulation studies. The first two studies have fixed 
bending mode frequencies centered at 6.0 and 10.0 hertz. 

The third study has the bending mode frequencies again fixed 
but at 7.0 and 9.0 hertz to ensure that the algorithm can 
separate two close bending mode frequencies. The fourth 
study has the bending modes fixed at 5.5 and 14.0 hertz to 
test the ability of the algorithm to properly filter two 
widely separated bending mode frequencies. The final two 
studies have dynamic bending modes with the upper bending 
mode changing at the rate of 1.5 hertz per second and the 
lower bending mode changing at 1.0 hertz per second. The 
rates of change far exceed the rates of change anticipated 
but were included as a demonstration of the capability of 
the algorithm to track very high rates of change in the 
bending mode frequencies. The final two studies illustrate 
the response of the system with the bending mode frequencies 
converging and diverging. 

Table 4-1 provides a summary of the conditions used in 
the six simulation studies. 

D. SIMULATION RESULTS 

Figures 4-2 and 4-5 illustrate the system response for a 
0.1 radian step pitch command with the bending mode fre- 
quencies fixed at 6.0 and 10.0 Hz. Line 1 in Figure 4-2 is 
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TABLE 4-1 



FIGURES 


INPUT 


BENDING MODES 


4-2 


0 . 1 rad 


6.0 Hz and 10.0 Hz 


4-3 


STEP 


FIXED 


4-4 


0.05*T rad/sec 


6.0 Hz and 10.0 Hz 


4-5 


RA.AIP 


FIXED 


4-6 


0.1 rad 


7.0 Hz and 9.0 Hz 


4-7 


STEP 


FIXED 


4-8 


0.1 rad 


5.5 Hz and 14.0 Hz 


4-9 


STEP 


FIXED 


4-10 


0.1 rad 


6.0 Hz + 1.0 Hz/sec 


4-11 


STEP 


DYNAMIC 


4-12 


0 . 1 rad 


8.0 Hz - 1.0 Hz/sec 


4-13 


STEP 


9.0 Hz + 1.5 Hz/sec 






DYNAMIC 
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Time, Sec. 



System; Third order missile model with fixed 
bending modes at 6 and 10 Hz 

Input; 0.1 radian step command 

Line 1 is the nominal missile response without 
bending modes 

Line 2 is the missile response with bending modes 

Figure 4-2 Missile Step Response with Fixed 6 and 
10 Hz Bending Modes 
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Figure 4-3 Frequency Detector Output for Fixed 
6 and 10 Hz Bending Modes 
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the nominal response of the plant without bending modes to 
the same input. It is provided as a standard by which the 
effectiveness of the adaptive digital notch filtering may 
be evaluated. Line 2 in Figure 4-2 is the response of the 
system with the bending modes and the adaptive digital notch 
filters , 

Figure 4-3 provides a history of the center frequency of 
the high and low channel notch filters. Line 1 is the his- 
tory of the low channel C5. 0-8.0 Hz) notch center frequency 
and Line 2 is the history of the high channel (8.0-12.0 Hz) 
notch center frequency. Lines 3 and 4 indicate the actual 
bending mode center frequencies. 

Figure 4-2 must be compared with Figure 2-17 to gain the 
proper perspective on the degree to which the system has been 
stabilized by the adaptive notch filters. Figure 2-17 is 
precisely the same system and the same input except the 

I 

adaptive digital notch filters were not in the system simu- 
lated in Figure 2-17. The amount of leakage of the 6.0 
hertz bending mode frequency through the notch filer is due 
to the relatively low precision coefficients used in the 
notch filters. Simulation with double precision coefficients 
reduced the leakage considerably but the implementation of 
the notch filter with this precision was not considered 
practical for implementation with a microprocessor and 
therefore not pursued. 
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Figure 4-3 illustrates that the algorithm can rapidly and 
accurately track the 6.0 hertz bending mode frequency but 
exhibits limit cycle behavior when tracking the 10.0 hertz 
bending mode. The limit cycle behavior can be attributed to 
the low frequencies associated with the plant dynamics causing 
errors in the frequency detector as discussed in Chapter 
Three. The limit cycle does not appear to diminish the high 
channel's effectiveness in eliminating the 10.0 Hz bending 
mode in that there is no evidence of the 10.0 hertz bending 
mode in the output and it is attenuated to the point that the 
frequency detector can no longer track it after 1.3 seconds. 

Figures 4-4 and 4-5 show the response of the system for 
a 0.05*T ramp pitch command with bending modes centered at 
6.0 and 10.0 hertz. Comparing Figure 4-4 with Figure 2-23 
it is evident that when the adaptive digital notch filters 
are included in the system the response is essentially the 
same as if no bending modes were present. Figure 4-5 shows 
the ability of the high channel to much more accurately lock 
onto the 10.0 hertz bending mode when the high plant rates 
are not present. The very satisfactory performance for a 
ramp input indicated that this input did not stress the 
ability of the adaptive digital notch filter algorithm and 
it was therefore not used in the remainder of the simulation 
studies . 

Figures 4-6 and 4-7 show the response of the system for 
a 0.1 radian step input when the bending mode frequencies 
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System; Third order missile model with fixed 
bending modes at 6 and 10 Hz 

Input; 0.05 radian/sec ramp command 

Line 1 is the nominal missile response without 
bending modes 

Line 2 is the missile response with bending modes 

Figure 4-4 Missile Ramp Response with Fixed 6 and 10 Hz 
Bending Modes 
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Figure 4-5 Frequency Detector Output for Fixed 
6 and 10 Hz Bending Modes 
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System: Third order missile model with fixed 

bending modes at 7 and 9 Hz 

Input: 0.1 radian step command 

Line 1 is the nominal missile response without 
bending modes 

Line 2 is the missile response with bending modes 



Figure 4-6 Missile Step Response for Fixed 7 and 
9 Hz Bending Modes 
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Figure 4*7 Frequency Detector Output for Fixed 
7 and 9 Hz Bending Modes 
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were fixed and only separated by 2.0 hertz. The system 
response is very close to the nominal plant with some evi- 
dence of both bending modes present in the output. The 
cross channel leakage can be seen to cause oscillations in 
both of the frequency detectors. 

Figures 4-8 and 4-9 show the response of the system to 
a 0.1 radian step input ivith the bending modes separated by 
7.5 hertz. The lower frequency bending mode is at 5.5 hertz 
or 0.5 hertz above the lower limit to which the algorithm 
is alloi\fed to track and the upper frequency bending mode is 
14.0 hertz. The upper bending mode frequency is outside of 
the passband for the high channel frequency detector but 
within the 15.0 hertz that the high channel is allowed to 
track a bending mode. The performance of the system is 
marginal at best. This simulation clearly illustrates the 
inability of the algorithm to reliably acquire or track a 
bending mode higher than 12.0 hertz. The problem is not 
accentuated by the bandpass filter in the high channel 
attenuating frequencies above 12.0 hertz. Increasing the 
bandwidth of the high channel to 15.0 hertz did not enhance 
the ability of the algorithm to reliably identify a fre- 
quency above 12.0 hertz. 

Figures 4-10 and 4-11 show the system response for a 0.1 
radian step input with dynamic, converging bending mode 
frequencies. The high frequency bending mode starts at 12.0 
hertz and decreases at 1.5 hertz per second to 9.0 hertz. 
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System: Third order missile model with fixed 

bending modes at 5.5 and 14.0 Hz 

Input: 0.1 radian step command 

Line 1 is the nominal missile response without 
bending modes 

Line 2 is the missile response with bending modes 



Figure 4-8 Missile Response with Fixed 5.5 and 14.0 Hz 
Bending Modes 



113 



Q 

O 



Frequency 

Hz 




Figure 4-9 Frequency Detector Output for Fixed 
5.5 and 14.0 Hz Bending Modes 
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System; Third order missile mode with converging 
dynamic bending modes 

Input: 0.1 radian step command 

Line 1 is the nominal missile response without 
bending modes 

Line 2 is the missile response with bending modes 



Figure 4-10 Missile Step Response with Converging 
Dynamic Bending Modes 
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Figure 4-11 Frequency Detector Output for Converging 
Dynamic Bending Modes 
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The low frequency bending mode starts at 6.0 hertz and in- 
creases at 1.0 hertz per second to 8.0 hertz. The ability 
of the adaptive digital notch filters to stabilize the system 
with such high rates of change was not expected. When 
higher rates of change were simulated for the converging 
case, the error in the high frequency estimate became so poor 
that the system became unstable at the frequency of the high 
bending mode. A rate of 1.5 hertz per second is felt to be 
the limit on the rate of change of either bending mode fre- 
quency. 

Figures 4-12 and 4-13 show the system response to a 0.1 
radian step input with dynamic diverging bending mode fre- 
quencies. The upper bending mode frequency starts at 9.0 
hertz and increases at the rate of 1.5 hertz per second to 

12.0 hertz. The lower frequency bending mode starts at 8.0 
hertz and decreases at the rate of 1.0 hertz per second to 

6.0 hertz. The very good performance of the system is pro- 
bably due to the notch filters being initialized at 8.0 
hertz. The inability of the high channel to track a bending 
mode frequency greater than 12.0 hertz is evident as the 
high bending mode frequency increases. 
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System: Third order missile model with diverging 

dynamic bending modes 

Input: 0,1 radian step command 

Line 1 is the nominal missile response without 
bending modes 

Line 2 is the missile response with bending modes 

Figure 4'-12 Missile Step Response with Diverging 
Dynamic Bending Modes 
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Figure 4-13 Frequency Detector Output for Diverging 
Dynamic Bending Modes 
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V. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

The feasibility of implementing an adaptive digital notch 
filter to stabilize the Trident Missile Control System when 
two dynamic bending modes are present has been demonstrated. 
The conditions simulated in Chapter Four represent the worst 
case conditions and the performance of the algorithm under 
these conditions was satisfactory considering the very de- 
sirable response to less violent and more likely commands 
as shoivn in Figure 4-4. 

The algorithms used in the simulation studies contained 
in Chapter Four were developed with the ability to implement 
them in a microprocessor program as a primary concern. Sim- 
plicity of the algorithm was maintained, perhaps at the 
expense of performance, to ensure that the algorithm could 
be implemented on a microprocessor for subsequent real-time 
studies. The time varying, nonlinear characteristics of the 
system when the adaptive digital notch filters are tracking 
dynamic bending modes with varying input commands deny an 
analytical approach and therefore require extensive simula- 
tion to optimize the algorithm. 
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B. RECOMMENDATIONS FOR FURTHER STUDIES 



The simplicity o£ the digital bandpass filters allows an 
excessive amount of out of band interference. The perfor- 
mance of the bandpass filters could be significantly improved 
without an appreciable increase in the order of the filters 
by using more sophisticated types of filters such as a 
Chebyshev or an elliptical filter. The temptation to enhance 
the performance of the bandpass filter in nonreal-time simu- 
lations without considering the implementation problem on a 
microprocessor must be resisted. 

The direct realization of the bandpass and notch digital 
filters did not minimize memory, noise due to finite pre- 
cision coefficients or the number of multiplications required 
for each output. The structure of the filters needs to be 
examined with a specific microprocessor's capabilities con- 
sidered. The ZILOG Z-8000 microprocessor capabilities were 
considered during the development of the algorithm as a cur- 
rent state-of-the-art microprocessor. 

The algorithm must be evaluated on the six degree of 
freedom model of the Trident Missile available at Lockheed 
Missile and Space Company, Sunnyvale, California. The 
results of the simulation studies contained in Chapter Four 
are only as valid as the assumptions made in modeling the 
bending modes and in simplifying the plant. It must be 
emphasized that the algorithm contained in this thesis will 
not perform satisfactorily with the bending mode frequencies 
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present in the current Trident Missile. The algorithm in 
this thesis was designed to identify, track and attenuate 
bending modes that are a full octave below those of the 
current Trident Missile. 

A phaselock loop scheme should be reconsidered as a 
frequency detector. The phaselock loop has precisely the 
desired characteristic of seeking out and locking onto a 
frequency within its predetermined acquisition range and 
seems ideally suited for the task of identifying and track- 
ing the bending mode frequencies. It is recommended that 
future studies of the application of phaselock loops in 
this application use actual phaselock loops rather than com- 
puter simulation due to the inability to properly simulate 
the non linear behavior of the phaselock loop with algorithms 
known to the author. The feasibility of prefiltering the 
input to the phaselock loop with a bandpass filter to atten- 
uate out of track range signals and thereby reduce the lock- 
on time should be considered. 
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INITIAL DIGITAL BENDING MODE FILTER PARAMETERS 
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HZ DIGITAL BANDPASS FILTER AND FREQUENCY DETECTOR 
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COMPUTE THE CENTER FREQUENCY FOR THE HIGH FREQUENCY NOTCH 
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MODEL OF THE THIRD ORDER PLANT 
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ASSOCIATED FIGURES 
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INITIAL DIGITAL BENDING MODE FILTER PARAMETERS 
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HZ DIGITAL BANDPASS FILTER AND FREQUENCY DETECTOR 
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COMPUTE THE CENTER FREQUENCY FOR THE HIGH FREQUENCY NOTCH 
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HZ DIGITAL BANDPASS FILTER AND FREQUENCY DETECTOR 
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COMPUTE THE CENTER FREQUENCY FOR THE HIGH FREQUENCY NOTCH 
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MODEL OF THE THIRD ORDER PLANT 
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THET 



INITIAL DIGITAL BENDING MODE FILTER PARAMETERS 
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HZ DIGITAL BANDPASS FILTER AND FREQUENCY DETECTOR 



OsJ 

CM 



CD 

o 

vO 

00 

CO 

• 

0 

1 



<NJ 

CD 

(M 

OD>t‘ 

CNJ CD 

t < 

^00 

+ • 

(\iO 

^ I 



cD>r 

incD 

mm 

0 (\l 

mr^ 
O'i- 
• vO 

O t 



-jm 

cD^ 

^CD 

m • 

(NJO 

^ I 
•h- 
03 

I a 

CD^ 

•tt'CD 

m-J> 

• 



3 3 3 

eg ^ CO i-< si" CNJ 



OO*- 
l) It 

00 HKCDCOCDCDCDCDCDCO 

33 II II II II II II II It 
O OOrg.^<\l*^<M^(\j^ 

^fNi.^.- 4 <M<MmmNt<f 



m 



COCDCDCDCQCQCDCDCQCD 



CNJ 

CM 

X 

CO 

m 



Qc: 


fM 


• 


o 


CNJ 


X 


♦- 


h- 


o 


o 


• 


h- 


LU 


o 


o 


1 - 


1 


2 


LU 






Q 


CM 


> 




XfM 


O 


V 


CD>t 


2 


O 




LU 


2 


mcD 


3 


LU 


ms^ 


a 


3 




LU 


a 


moo 


Qc: 


LU 


• • 


LL 


a: 


-40 




u. 


-f 1 


2 




CM^ 


O 


Q 






2 


XX 




< 


CD CD 


LU 






X 


a: 


(Min 


1 - 


LU 


vOrsj 




h— 




CC 






o 






LL 


LL 


o • 








> 


00 


—4 + 


O 


C/^ 


^(M 


2 


< 


xm 


UJ 


Q. 


CDX 


3 


Q 


* 1 ^ CO 


Of 


2 




LU 


< 


vQ ^4 


a: 


CD 


• 


LL 




mo 






• 1 


Qi: 


< 


o>- 


LU 


h- 


1 3 


1 - 




DlO 


2 


O 


rQ 


LU 




^l*X 


O 


Q 


in CD 






vO-J^ H- H- 1 - 


LU 


rsj 


3 3 3 


X 


X 


• • ^ *-4 c 1 — » O •— 40 


K 




00-4 <Mr 4 m- 4 NtCM 




m 


II II xcx:xxxxxx 


LU 


^4 


h- H- CD CD CD CD CD CD CD CD 






33 II II It II II II II II 


3 


o 


O O CM CM CM -4 CM -4 


CL 


H* 


^ CM ^ ^ rgcM cn m ' 4 " 


2 




xxxxxxxxxx 


O 


CO 


CD CO CO CD CD CO CD CD CD CQ 


o 



(>0 

3 

Q 

< 

0 ^ 







LU 






Q 






•- 




o 


a 




m 


LU 






QC 


o 


oo 


LL 


• 


mh- 


w 


oo 




IL 


•iin 


OO 


LL 


h- 


KO O 


LU 


30 


— »• • 


O 


Ol- 


0—4 CM 


O 


1-4 




II 


•JO 


CMZ'- 


«-4 


CDO*^ 


*-^v 0 »-^O 


2 


1 — 


^ 0 ^-*UJ 


N -4 


LU— 


• •■«• d: 


< 


XO »-4 


000 LL .-40 


1-4 •»— 


• • • + OO ^ 


l-O-J 


l-h-(Ml~ 3 <M 



w • I --lO-jwUJQ^a 



000^0 • tV..</)<CDUJ 

ooujoK^^o^Qc: •‘o: 

m •!- II •11- II ^LL 

CC^ II It 

II II 3CDH 
II h-KZH-H»C3ar»-t P^XI 
,-4 LU O ^00 

O LLU-QC QC <5 -4 

K I-Il- U •-•U- u. oc: < LL 



* -IV 






•». -II. - 1 ^ 






202 



COMPUTE THE CENTER FREQUENCY FOR THE HIGH FREQUENCY NOTCH 
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MODEL OF THE THIRD ORDER PLANT 
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EL OF THE PLANT WITHOUT BENDING MODES 
ESTABLISH NOMINAL PERFORMANCE . 
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ASSOCIATED FIGURES 
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INITIAL DIGITAL BENDING MQDE FILTER PARAMETERS 
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HZ DIGITAL BANDPASS FILTER AND FREQUENCY DETECTOR 
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COMPUTE THE CENTER FREQUENCY FOR THE HIGH FREQUENCY NOTCH 
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MODEL OF THE THIRD ORDER PLANT 
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INITIAL DIGITAL BENDING MODE FILTER PARAMETERS 
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COMPUTE THE CENTER FREQUENCY FOR THE HIGH FREQUENCY NOTCH 
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MODEL OF THE THIRD ORDER PLANT 
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